Qiu et al. BMC Genomics 2012, 13:397 
http://www.bionnedcentral.conn/1471 -21 64/1 3/397 



Genomics 



RESEARCH ARTICLE Open Access 



Comparative transcript profiling of gene 
expression between seedless Ponkan mandarin 
and its seedy wild type during floral organ 
development by suppression subtractive 
hybridization and cDNA microarray 

Wen-Ming Qiu, An-Dan Zhu, Yao Wang, Li-Jun Chai, Xiao-Xia Ge, Xiu-Xin Deng and Wen-Wu Guo" 



Abstract 

Background: Seedlessness is an important agronomic trait for citrus, and male sterility (MS) is one main cause of 
seedless citrus fruit. However, the molecular mechanism of citrus seedlessness remained not well explored. 

Results: An integrative strategy combining suppression subtractive hybridization (SSH) library with cDNA microarray 
was employed to study the underlying mechanism of seedlessness of a Ponkan mandarin seedless mutant {Citrus 
reticulata Blanco). Screening with custom microarray, a total of 279 differentially expressed clones were identified, 
and 133 unigenes (43 contigs and 90 singletons) were obtained after sequencing. Gene Ontology (GO) distribution 
based on biological process suggested that the majority of differential genes are involved in metabolic process and 
respond to stimulus and regulation of biology process; based on molecular function they function as DNA/RNA 
binding or have catalytic activity and oxidoreductase activity. A gene encoding male sterility-like protein was highly 
up-regulated in the seedless mutant compared with the wild type, while several transcription factors (TPs) such as 
AP2/PRPBP, MYB, WRKY, NAG and G2G2-GATA zinc-finger domain TPs were down-regulated. 

Conclusion: Our research highlighted some candidate pathways that participated in the citrus male gametophyte 
development and could be beneficial for seedless citrus breeding in the future. 
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Background 

Seedlessness is a desired fruit trait for consumers, and a 
fruit is considered to be seedless if it produces no seeds, 
traces of abortion seeds, or significant reduced-number of 
seeds [1]. Some plants can set seeds asexually through 
apomixis. However, in most flowering plants, seed initi- 
ation requires signals activated by the double fertilization 
event that occurs in the embryo sac, and seed is produced 
sexually from the fertilized ovule [2,3]. Various phytohor- 
mones such as gibber ellins (GAs), auxins and cytokinins 
are involved in this signaling process [4-6]. GAs and 



* Correspondence: guoww(5)mail. hzau.edu.cn 

Key Laboratory of Horticultural Plant Biology (Ministry of Education); National 
Key Laboratory of Crop Genetic Improvement, Huazhong Agricultural 
University, Wuhan 430070, China 

(3 BioMed Central 



jasmonic acid/jasmonate derivatives (JAs) were found to 
play crucial roles in plant reproductive development [7,8]. 

Citrus is one of the most important fruit crops with great 
economic and health value around the world [9]. However, 
some citrus varieties are seedy, and seedy fruits have con- 
strained the development of fresh citrus market. Therefore, 
breeding seedless citrus varieties is a long-term pursuit for 
citrus breeders worldwide [10,11]. Nowadays, Satsuma 
mandarin and navel orange are two of the most famous 
and widely grown citrus varieties, mainly due to their seed- 
less trait. For decades, great progress on seedless citrus 
breeding was made by traditional approaches such as sex- 
ual hybridization, seedling and bud sport mutation. How- 
ever, due to the peculiarities of citrus reproductive biology 
such as long juvenile period and nucellar polyembryony. 
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traditional breeding is inefficient and costly [12]. Modern 
biotechnological approaches (e.g. somatic hybridization) 
have potential to effectively expedite breeding process of 
citrus [13-15]. As most citrus varieties can produce fruits 
parthenocarpically [16], male or female sterility, embryo sac 
abortion, self-incompatibility, polyploidy and even environ- 
mental stress can result in seedless citrus fruits [17,18]. Ac- 
tually there were some successful reports about seedless 
fruit production by genetic transformation. Ectopic expres- 
sion of iaaHgene with DefH9 as promoter to elevate auxin 
levels in placenta or ovules resulted in seedless fruits 
[19,20]. Another effective strategy was by specific expres- 
sion of toxin proteins during early development of plant 
reproductive organs. Typical cases were the ectopic trans- 
formation of the Barnase gene from Bacillus amyloliquefa- 
ciens [21,22]. Potential cases were by specific expression of 
enzymes such as chloroplast Chaperonin 21 and ubiquitin 
extension protein S27a to induce cell disruption of seed 
tissues for parthenocarpic plants [11,23,24]. And in our 
laboratory, the Arabidopsis thaliana MAC12.2 gene had 
been introduced into precocious trifoliate orange {Poncirus 
trifoliata [L.] Raf) for production of potential seedless 
fruits [25]. 

Male sterility (MS) is one of the main causes for seed- 
less fruit production in citrus. In recent years, great pro- 
gress on MS was made with annual plants especially rice 
[26,27], Arabidopsis [28] and oil-rape [29], and a serial 
of genes regulated tapetum, anther and pollen develop- 
ment were identified. However, there remained very 
limited information on MS of perennial woody plants 
such as citrus. Ponkan mandarin {Citrus reticulata 
Blanco) is a widely grown citrus variety in China. Within 
this variety, many variants were derived through sexual 
hybridization and mutation such as bud sport mutation. 
'Qianyang seedless' Ponkan mandarin (QS) is an elite 
seedless variant selected from bud sport mutation of a 
common seedy Ponkan mandarin, and it can set fruits 
with no seeds (even no seed rudiments) in open orchard 
[30,31]. In this article, QS and a common seedy Ponkan 
mandarin 'Egan NO.l' (EG) were used for comparative 
study. These two mandarins shared highly close genetic 
relationship based on molecular marker analysis and 
showed no distinctly morphological differences except 
that QS was completely male sterile while Egan No 1 has 
normal flower. In order to gain general understanding 
on genes involved in this MS mutation, suppression 
subtractive hybridization (SSH) [32] combining with 
cDNA microarray was performed to detect differen- 
tially expressed genes. Several candidate genes and 
related pathways were focused in particular. Our re- 
search identified some useful genes which could be 
beneficial to citrus seedless breeding. The results 
could help to reveal the molecular mechanism of 
male sterility of Ponkan mandarin and shed light on 



seedless trait formation of other perennial woody plant at 
the gene expression level. 

Results 

Phenotype analysis of the floral organs of QS 

Previous studies suggested that the floral organs (actu- 
ally the whole plant) of QS had no morphological diffe- 
rence from the wild type. To further validate the 
phenotype of this seedless Ponkan mandarin, we mea- 
sured the length of filament and pistil, and the average 
ratio of filament to pistil (filament length/pistil length) 
was 0.83 ± 0.01 for EG and 0.79 ± 0.01 for QS. And for 
EG, the pistil was 0.155 ± 0.01 cm longer than filament 
while for QS, the pistil was 0.166 ± 0.009 cm longer than 
filament. Above data further confirmed that the floral 
organs of both EG and QS had no morphological differ- 
ence, and the seedless trait was not caused by malforma- 
tion of reproductive organs. However, the number of 
pollen grains per anther of QS was 9.5% less than that of 
EG. The pollen dying viability of QS was 6.0% ± 1.0% 
(or 6.5% ± 1.0% for I2-KI2 staining) in striking contrast 
to the high viability of 93.8% ± 0.9% (or 89.6% ± 2.5% 
for I2-KI2 staining) for EG. Pollen germination test 
found that no pollen of QS could germinate. Further- 
more, SEM assays showed abnormal structures of the 
pollen grains of QS (Figure 1), confirming that QS is 
male sterile. 

Construction of SSH-cDNA libraries and overall feature of 
the differential transcript profiling 

To identify genes associated with the MS of QS, SSH 
cDNA libraries (both forward and reverse) were con- 
structed from floral organs of QS and EG. A total of 
6,048 cDNA clones derived from the SSH-cDNA librar- 
ies including 4,195 from the forward library and 1,853 
from the reverse one were successfully amplified, and 
then used for a custom cDNA microarray. Each cDNA 
clone has triplicate spots on the array. The RNA samples 
of the four developmental stages (SF, MF, BF and OV) 
were used for array-hybridization. The fluorescent dye- 
labelled cDNA and hybridization strategy was employed 
for the microarray assay. 

From the 6,048 clones printed on the glass slide, 279 
cDNA clones (278 non-redundant) were differentially 
expressed (false discovery rate (FDR) <0.05 and a fold 
change > 2) between QS and EG. Among these cDNA 
clones, 218 (78%) were down-regulated while only 61 
(22%) showed up-regulated expression across the four de- 
velopmental stages; and the differentially expressed clones 
peaked at full bloom stage (BF) (Figure 2). At this stage, 
many more clones showed down- regulated than up- 
regulated expression. During the four developmental 
stages, one clone (GenBank accession no. JU497336) en- 
coding a putative cysteine protease (tr[B9RRA4]) showed 
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Figure 2 Number of clones significantly up- and down-regulated in QS during four developmental stages. Each two numbers in 
parentheses of x-axis indicated the numbers of clone up- and down-regulated respectively. 



down-regulated expression at BF stage but up-regulated at 
OV stage (young ovaries of 2-3 days after flowering). 

Sequencing of the differentially expressed clones and EST 
analysis 

Among the 279 difl'erentially expressed clones, 255 non- 
redundant clones were subjected to one single-pass se- 
quencing. In all, 237 high-quality ESTs (average length 
was 496 bp) were yielded after eliminating vectors and 
unreliable sequences. These ESTs were assembled using 
CAP3 program, and 133 unigenes (43 contigs and 90 
singletons) were obtained with sequence redundancy of 
43.9%. The majority of the contigs (38) contained 2-5 
ESTs, whereas only 5 contigs contained 6-11 ESTs, indi- 
cating an ideal normalization and subtraction. Of the 
133 unigenes, 80 (60.1%) showed differential expression 
at BF stage. Subsequently, BLASTX search of the UniProt 
database showed that 20 unigens (15.0%) did not have sig- 
nificant hits (E-value <1.0 x e"^). However, when the 20 
unigenes were used in BLASTN (E-value <1.0 x e"^^) 
search of the Citrus Clementina transcript database [33,34] 
with local Blast software (ftp://ftp.ncbi.nlm.nih.gov/blast/ 
executables/release/LATEST/), 17 genes had significant 
hits and high scoring pairs (HSP) showed high nucleotide 
identity. It suggested that these 20 unigenes were unique 
for citrus, and three of them were novel citrus genes. 

Based on the microarray analysis, the relative expres- 
sion profiles of all 255 ESTs were performed hierarchical 
clustering with cluster software (version 3.0). Four typ- 
ical relative expression patterns were observed in QS 
versus EG at four developmental stages. Figure 3 A and 
3B showed a group of clones down-regulated mainly at 



squaring stage (SF) and full bloom stage (BF), respect- 
ively, while the other two groups of clones were down/ 
up-regulated constitutively during the developmental 
stages (Figure 3C and 3D). In addition, candidate genes 
with putative function that could be important for the 
MS of QS were specifically collected (Table 1). It is note- 
worthy that 27.7% of the unigenes (not listed in the 
table) were only annotated as putative proteins or with 
no defined biological process besides 15% unigenes with 
no hits in the database. 

GO annotations were conducted and three categories 
representing molecular functions, biological processes, 
and cellular components were assigned. Figure 4 showed 
the percentage distributions of GO terms (2nd level GO 
terms) based on biological process. It indicated that dur- 
ing the floral organ development, the majority of differ- 
entially expressed genes were involved in metabolic 
process (46%) or responded to stimulus (27%) and regu- 
lation of biological process (18%). In addition, the other 
two GO categories (molecular functions and cellular 
components) were also generated (data not shown). In 
the molecular function category, large proportion of uni- 
genes may have binding activity (59%), catalytic activity 
(19%), or oxidoreductase activity (11%), while the cellu- 
lar components consisted mainly of intracellular (57%) 
and membrane (23%). 

Metabolic pathways involved in formation of seedless 
fruit 

As large proportion of altered expressed genes were 
involved in varieties of metaboUc processes. Based on 
the KEGG (Kyoto Encyclopedia of Genes and Genomes) 
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Figure 3 (See legend on next page.) 
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(See figure on previous page.) 

Figure 3 Cluster analysis of expression profiles of altered expressed genes in the QS versus EG. A sliowed a cluster of ESTs tliat were 
down-regulated mainly at squaring stage (SF) when the tetrads were produced and the microsporocyte underwent meiosis. B showed ESTs that 
were down-regulated especially at full bloom stage (BP). C and D suggested a cluster of ESTs that were down-regulated and up-regulated 
constitutively during the four developmental stages respectively. The ratio value was log2 transformed for each gene and used for the 
hierarchical clustering analysis. 

V J 



analysis, 36 different metabolic pathways were altered 
during the four developmental stages. Among these 
pathways, nine (25%) were related to amino acid meta- 
bolic pathway (Table 2), and genes involved in carbohy- 
drate and energy metabolism showed down- regulated 
expression during subsequent developmental stages of 
floral organs. Besides, genes related to specific secondary 
metabolism such as terpenoids and polyketides metabo- 
lism were also found to be altered. Interestingly, a gene 
(JU497309) encoding fatty acyl-CoA reductase, which 
may be involved in lipid metabolic process, was identi- 
fied (Table 1). This gene was found highly homologous 
with putative male sterile protein (GI: 255576327) in 
castor bean, fatty acyl-CoA reductase 3 (GI: 359500474) 
in poplar and male sterile 2-like protein {MS2) (GI: 
3549681) in Arabidopsis. Herein, this gene was named 
as male sterile-like protein. And qRT-PCR analysis 
showed its expression level increased from SF to BF 
stages and then declined at OV stage. The expression 
pattern was similar in both QS and EG; however, it 
showed obviously higher expression level in QS than in 
EG during the developmental process (Figure 5). 

Differential expression of transcription factor genes 

It is noteworthy that among the 133 unigenes, 12 were 
assigned to the category of transcription factor (TF) 
based on plant TF database (http://planttfdb.cbi.edu.cn/). 
Figure 6 showed the specific expression pattern of six 
AP2-ERF family TFs, two zinc-finger TFs, one MYB TF 
and one NAG TF using qRT-PCR assay. These TFs (ex- 
cept of NAC TF) had similar expression profile during the 
four developmental stages between EG and QS. For in- 
stance, among six AP2-ERF TFs, four (AP2-EREBP TFl, 
AP2-EREBP TF3, AP2/ERF domain containing TF2 and 
CBF/DREB-like TF) showed co-expression pattern like 
"V" type. It showed that the gene expression level in QS 
was higher than that in EG from SF stage to MF stage; 
however, these genes were subsequently repressed more 
obviously in QS from MF stage to BF stage, and the gene 
expression level was down-regulated mostly at BF stage. 
Two zinc-finger TFs (GATA TFS and Cys2-His2 type) and 
one R2R3-MYB TF likewise showed similar "V" type- 
variation tendency. The other two AP2-ERF TFs (AP2- 
EREBP TF2 and AP2/ERF domain containing TFl) 
showed "V"-like type expression pattern in QS. However, 
the expression pattern of AP2/ERF domain containing 



TFl was somehow different from others, as it showed 
relatively stabilized expression level during the four stages 
in EG. As for NAC TF, its expression level was down- 
regulated obviously at BF and OV stages in QS compare 
with EG. It was notable that no expression was observed 
at OV stage in QS. The results suggested that these TFs 
could play important roles in the seedless phenotype for- 
mation, and the relative expression level in QS versus EG 
seemed to be key factor in this process. 

Verification of microarray data 

Two approaches were used to examine the quality of the 
microarray data. First, as one contig was assembled by 
several ESTs that were arrayed at random location in the 
microarray, so these ESTs sharing similar sequence or 
encoding the same gene would share similar expression 
pattern. Additional file 1: Figure SI showed that four 
ESTs (F2-13 G, F6-15I, F7-180, F8-12A) were assembled 
into one unigene (JU497321) which encoded methionine 
synthase, and these four ESTs truly shared similar ex- 
pression pattern. For the other approach, qRT-PCR was 
performed on 11 unigenes using gene-specific primer 
pairs. Expression patterns were compared at the four 
developmental stages between QS and EG. Additional 
file 2: Figure S2 showed the correlation analysis of the 
ratio values of differential expression level from micro- 
array to that from qRT-PCR. Linear regression [(average 
microarray ratio value) = a (qRT-PCR value) + b] analysis 
showed a good coefficient of variation {R^ = 0.847). 
These results confirmed the reliability of the microarray 
data. 

Discussion 

Here, we combined SSH and microarray techniques to 
investigate potential mechanism underlying seedlessness 
in Ponkan mandarin. SSH was proved to be an efficient 
and popular approach to enrich and identify differen- 
tially expressed genes between wild-type and its mutant 
or treatment [35,36]. However, because of high sensitiv- 
ity of SSH, usually a large number of clones could be 
obtained but inevitably included some false-positive 
ones. Screening the SSH libraries to identify some candi- 
date genes using microarray and to validate using qRT- 
PCR has proved to be a high-throughput and efficient 
way [37-39]. However, relatively few clones were isolated 
in this study. Of the 6,000 clones, only 279 cDNA clones 
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Table 1 List of selected candidate functional genes related to the formation of the phenotype of QS 

GenBank accsession no. EMI ^ Description ^ e-Value ^ Clones ^ 

Up-regulated 



JU497309 


B9SST7 


Male sterility protein 


2E-77 


7 


JU497311 


B9S7V4 


STS14 protein 


6E-23 


3 


JU497315 


Q5CD81 


(E)-beta-ocimene synthase 


lE-120 


2 


JU497324 


B9RXQ0 


Tryptophan synthase beta chain 


6E-91 


2 


JU497327 


B9RT81 


Multidrug resistance pump 


5E-50 


3 


JU497333 


B9S7Q2 


Zinc finger protein 


lE-67 


3 


JU497348 


A9ZN18 


Geranyl-diphosphate synthase 


3E-11 


5 


JU497417 


B9R726 


Xylem serine proteinase 1 


7E-24 


1 


JU497359 


B9SCQ9 


Oxidoreductase 


2E-36 


1 


JU497418 


A1ECJ7 


Putative miraculin-like protein 2 


7E-24 


1 


JU497422 


A9XCN2 


Putative DNA binding protein 


2E-24 


1 


JU497389 


B9I523 


Multidrug resistance protein ABC transporter family 


lE-121 


2 


JU497397 


B9SIR0 


S-locus-specific glycoprotein 


2E-45 


1 


JU497403 


B9S0U2 


Stomatin-1, putative 


lE-07 


1 


Down-regulated 










JU497308 


B9S7Q1 


Zinc finger protein 


lE-140 


2 


JU497318 


B9S1E9 


Transcription factor AtMYC2 


4E-20 


2 


JU497321 


B2VQE0 


Methionine synthase 


lE-107 


4 


JU497323 


P42802 


lnositol-3-phosphate synthase 


lE-112 


4 


JU497331 


B9SR02 


Multicopper oxidase 


6E-78 


2 


JU497332 


B3RGD7 


CBF/DREB-like transcription factor 


6E-44 


4 


JU497336 


B9RRA4 


Cysteine protease 


2E-25 


2 


JU497338 


Q8VWL8 


Beta-mannosidase 


6E-64 


2 


JU497342 


082547 


Chitinase CHIl 


8E-28 


7 


JU497343 


049817 


Late embryogenesis abundant protein 


2E-22 


2 


JU497344 


B9SI15 


Ubiquitin-protein ligase 


2E-54 


2 


JU497351 


A9PLA1 


AP2/ERF domain-containing transcription factor 


2E-46 




JU497352 


B9T724 


GATA transcription factor 


4E-62 




JU497353 


B9RZK6 


Protein COBRA 


4E-82 




JU497354 


A2IB54 


Mitogen-activated protein kinase 


4E-85 




JU497356 


B9RST2 


Glutamine synthetase plant 


3E-71 




JU497357 


Q3KN68 


Isoflavone reductase-like protein 5 


2E-56 




JU497361 


B9SQM6 


Transcription factor 


3E-34 




JU497362 


Q9FT21 


Putative glutathione S-transferase T3 


2E-07 




JU497364 


B9GWH5 


Glutathione peroxidase (Fragment) 


2E-25 




JU497367 


B9RN18 


Peptidyl-prolyl cis-trans isomerase 


3E-54 




JU497368 


B9R762 


R2R3-myb transcription factor 


lE-69 




JU497369 


082547 


Chitinase CHIl 


8E-28 




JU497419 


P83948 


Recti nesterase-3 


lE-46 




JU497372 


B9NBQ9 


AP2/ERF domain-containing transcription factor 


lE-108 




JU497373 


Q9ZRC9 


ACC oxidase 


6E-36 




JU497421 


A5YWA9 


NAC domain protein 


lE-51 




JU497378 


Q7Y066 


Plasma membrane H^-ATPase 


lE-114 





Qiu et al. BMC Genomics 2012, 13:397 
http://www.bionnedcentral.conn/1471 -21 64/1 3/397 



Page 8 of 1 7 



Table 1 List of selected candidate functional genes related to the formation of the phenotype of QS (Continued) 



JU497423 
JU497424 
JU497427 
JU497388 
JU497391 

JU497398 
JU497401 
JU497405 
JU497406 
JU497434 
JU497412 
JU497413 



B9RIP3 Hevamine-A 

B9T0Z2 ADP/ATP carrier protein 

Q84U94 Cliloropliyll a/b-binding protein (Fragment) 

A7XUL4 deliydration-responsive element binding protein 

Q4F6Y8 Putative AP2-binding protein 

Q6EV47 Non-specific lipid-transfer protein (Fragment) 

B9SJL5 Amine oxidase 

B9S925 Zinc finger protein 

B9T868 Putative peroxidase C3 (Fragment) 

Q8H2A1 CaffeoyI CoA 0-metliyltransferase (Fragment) 

B9HGW6 Glutaredoxin 

B9R762 R2R3-myb transcription factor 



lE-22 
lE-06 
2E-14 
lE-14 
6E-15 

lE-45 
lE-57 
3E-88 
8E-26 
3E-11 
6E-20 
3E-21 



^ The EMI codes of the most similar genes to the EST sequences. 
^ The description of sequences based on Uniprot database. 

The best e-value from a BLASTx search for corresponding EST sequences. 
^ The number of sequenced clones in the libraries. 

were identified as differentially expressed. Such results 
may suggest that there were little variations between QS 
and EG mandarins in gene expression. It was hypothe- 
sized that bud sport mutant was likely caused by single 



gene mutation, DNA methylation or retroelement activ- 
ity [40,41]. In this research, various types of DNA mar- 
kers including SCAR [42], and SSR (172 pairs of 
primers), MSAP (96 pairs of primers) and AFLP (13 



GO annotation for biological process 

0 10 20 



30 



40 



50 





cellular process 
multicellular organismal development 
nucleobase, nucleoside, nticleotlde and nucleic acid metatiolic process 

cell death 
cell communication 
regulation of biological process 
metabolic process 
cell differentiation - 
transport 

cellular amino acid and derivative metabolic process 
response to stimulus 



Percentage of ESTs 

Figure 4 Distribution of the unigenes according to the biological process (2nd level GO terms). A total of 86 unigenes were annotated As 
one gene product could be assigned to more than one GO term, the percentages will add up to more than 100%. 
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Table 2 List of differentially expressed genes involved in amino acid, carbohydrate, energy, terpenoid and polyketides 
metabolism based on KEGG pathway database 



KEGG 
pathways 



EC 

number 



Putative function 



QS/EG 
SF 



MF 



BF 



OV 



Amino acid metabolism 



JU497356 
JU497374 
JU497373 
JU497321 

JU497330 
JU497338 
JU497364 

JU497324 
JU497377 



6.3.1.2 
6.3.5.4 
1.14.17.4 
2.1.1.14 

4.2.1.78 
3.2.1.21 
1.11.1.1 

4.2.1.20 
2.4.1.12 



Carbohydrate metabolism 

JU497385 3.2.1.14 
JU497313 1.13.99.1 
JU497323 5.5.1.4 
JU497357 1.3.1.45 
Energy metabolism 
JU497406 1.11.1.7 
JU497378 3.6.3.6 
Terpenoids and polyketides 
JU497376 1.1.1.295 
JU497315 4.2.3.15 
JU497325 4.2.3.20 



Glutamate-ammonialigase 0.71 ±0.03 0.84±0.05 0.50±0.01 0.79±0.03 

Asparagine synthase 0.50±0.03 1.15±0.08 1.20±0.09 1.02±0.06 

Aminocyclopropane carboxylate oxidase 0.90±0.02 1.52±0.05 0.52±0.03 1.15±0.09 

5-nnethyltetrahydropteroyltriglutamate 1.62±0.1 1 0.97±0.08 0.50±0.02 0.90±0.08 
-homocysteine-methyltransferase 

(S)-norcoclaurine synthase 1.04±0.11 1.94±0.28 1.29±0.04 1.17±0.13 

Beta-glucosidase 0.94±0.08 0.47±0.04 0.18±0.01 0.88±0.17 

Phospholipid-hydroperoxide 0.53±0.04 0.97±0.10 0.97±0.05 1.08±0.01 
-glutathione peroxidas 

Tryptophan synthase 1.86±0.13 2.64±0.16 1.03±0.05 0.70±0.05 

lndole-3-acetatebeta-glucosyl transferase 0.51 ±0.01 1.51 ±0.04 1.25±0.02 0.48±0.04 

Chitinase 1.15±0.15 1.28±0.05 1.00±0.17 0.45±0.02 

Inositol oxygenase 1.04±0.09 1.08±0.04 0.79±0.07 0.49±0.03 

lnositol-3-phosphate synthase 1.11 ±0.08 1 .1 7±0.20 0.48±0.03 1 .42±0.1 3 

2{prime}-hydroxy isoflavone reductase 1 .1 0±0.05 1 .28±0.07 0.44±0.01 1 .1 7±0.04 

Peroxidase 0.99±0.05 1.02±0.01 1.12±0.09 0.54±0.06 

Proton-exporting ATPase 0.99±0.05 0.90±0.02 0.51 ±0.02 0.94±0.03 
metabolism 

Momilactone-Asynthase 1.42±0.54 1.09±0.20 1.05±0.03 1.10±0.10 

Myrcene synthase 1 .24±0.04 2.28±0.1 6 0.52±0.01 1 .1 0±0.1 1 

(R)-limonene synthase 1.27±0.18 0.91 ±0.03 0.50±0.02 1.10±0.17 



pairs of primers) were employed to analyze the poly- 
morphism between these two mandarins, and no repeat- 
able polymorphic bands were detected (data no shown). 
These results suggested that very few nuclear genes were 
altered during the developmental stages. 

For the four developmental stages we chose, immense 
efforts were taken to determine which time-point was 
pivotal for stamen development, but there has no criteria 
for citrus gametophyte development. Though criteria for 
gametophyte development was available in model plant 
Arabidopsis [43], it can not be directly applied herein. 
Semi-thin and paraffin sections were performed in this 
study to survey the microsporogenesis of QS, and it was 
found that abnormal tetrads produced at the tetrad stage 
and subsequently the microsporocyte underwent abnor- 
mal meiosis. This process mainly occurred at SF stage (the 
diameter of floral organs is about 3 mm) (unpublished 
data). Additionally, large proportion (about 59.7%) of dif- 
ferentially expressed genes was found in BF when the 



anthers and pollen grains were almost mature, indicating 
that this time-point might be also important. 

Amino acid metabolic process 

Of the metabolic pathways with altered expressed genes, 
25% were involved in amino acid metabolism. Amino 
acids were not only primary metabolic products for nor- 
mal growth and development but also cell signaling 
molecules and regulators of gene expression and protein 
phosphorylation cascade [44]. Interestingly, among these 
amino acid metabolism pathways, two genes were down- 
regulated across the developmental stages in QS versus 
EG, one (JU497356) encoding glutamate-ammonialigase 
(EC 6.3.1.2), the other (JU497338) encoding beta- 
glucosidase (EC 3.2.1.21). In higher plants, glutamate- 
ammonialigase catalyzes ATP-dependent conversion of 
glutamate and ammonia into glutamine which occupies 
a central position of amino acid metabolic pathway [45], 
and this metabolic process is critical for coordinating 
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Sample stages 

Figure 5 Relative expression (y-axis) of male sterile-like protein in QS versus EG during four developmental stages (x-axis) by qRT-PCR. 

Columns and bars represent the means and standard errors (n = 3) respectively. 



metabolic balance in rice [46]. And beta-glucosidase 
could be used for the cellulosic ethanol industry [47] 
and has diversity of functions in plants. In maize, Zm- 
p60J encoding a beta-glucosidase could release active 
cytokinin, and might function in vivo to supply the 
developing maize embryo [48]. Additionally, some beta- 
glucosidases affect the properties of cell wall [49] and 
are associated with freezing tolerance, such as the SFR2 
in Arabidopsis [50]. Some beta-glucosidases are related 
to the efficiency of microspore embryogenesis [51]. It is 
noteworthy that a gene (JU497374) encoding asparagine 
synthase (EC 6.3.5.4) was down- regulated exclusively at 
SF (early stage of stamen development). And asparagine 
is one central intermediate in nitrogen assimilation and 
transportation in plant [52,53]. Recent studies showed 
that this gene played important role in defense against 
pathogens and salt stress [54,55]. Additionally, genes 
related to carbohydrate metabolism and energy metabol- 
ism also showed down- regulated expression in QS 
mainly at BF and OV (late stage of stamen develop- 
ment). These results suggested that the vital activities of 
QS weakened during early development stages of sta- 
men, and the metabolic process of nutrition and energy 
was also impaired at subsequent stages of stamen devel- 
opment especially when the stamen was mature. 

Two genes involved in cysteine/methionine metabo- 
lism and participated in the biosynthesis of ethylene 
were also identified in this study. One (JU497321) encodes 
5-methyltetrahydropteroyltriglutamate-homocysteine S- 
methyltransferase (EC 2.1.1.14) is likely involved in the 
biosynthesis of L-methionine. And the methionine can be 
transformed into S-adenosylmethionine (SAM) (the 



precursor of ethylene) [56]. The other one (JU497373) 
encodes aminocyclopropane carboxylate oxidase (EC 
1.14.17.4) and is a pivotal enzyme during the biosynthesis 
of ethylene. In addition, genes involved in the synthesis of 
lAA (indole-3-acetic acid) were also identified such as a 
gene (JU497377) encoding Indole-3-acetatebeta-glucosyl- 
transferase (EC 2.4.1.121). These results implied that the 
endogenous phytohormones might be involved in the 
male gametophyte development of citrus. 

Transcription factors 

It was known that floral organ formation and function 
were influenced by TFs regulation. In our research, 
twelve unigenes were assigned to the category of tran- 
scription factor, and six of them were identified as AP2- 
ERF family members. AP2-ERF TF containing highly 
conserved AP2/ERF DNA-binding domain, is a large 
family unique in plant. In our research, four AP2-ERF 
members showed similar expression pattern. AP2- 
EREBP TFl was closely homologous with atERF107 
(AT1G19210). This gene was likely involved in the regula- 
tion of gene expression by stress factors and by compo- 
nents of stress signal transduction pathways. However, until 
now, no experimental evidence was available. AP2-EREBP 
TF3 showed high similarity with ERF5 (AT5G47230.1). 
ERF5 might play an important role in plant innate immu- 
nity likely through coordinating chitin and other defense 
pathways [57]. Other research suggested that ERF5 and 
ERF6 might potentially overlap in their function and 
acted as positive regulators of JA/ethylene-mediated 
defense [58]. In tomato, this gene was mainly involved in 
responses to drought and salt stresses [59]. As for AP2/ERF 
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Figure 6 (See legend on next page.) 
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(See figure on previous page.) 

Figure 6 Relative expression pattern of six AP2-ERF family TPs, two zinc-finger TPs, one MYB family TP and one NAC domain TP. Tine 
accession number of eacli TF was given inside tlie parentliesis. Relative expression was defined as tine expression level in QS versus EG. Columns 
and bars represent the means and standard errors (n = 3) respectively. 

V J 



domain containing TF2, its closest relative was ERF104 
(AT5G61600.1). Recent studies showed that ERF104 was 
in vivo substrate of MPK6, and ethylene could release 
ERF104 and allow liberated ERF104 to access target genes 
related to plant defense [60]. CBF/DREB-like TF was of 
high similarity with CBF4 (AT5G51990.1) which was crit- 
ical regulator involved in cold acclimation and drought 
adaptation [61,62]. 

In addition, AP2-EREBP TF2 was highly homologous 
with RAP2.4 (AT1G78080.1). RAP2.4 acted at or down- 
stream of a converging point of light and ethylene sig- 
naling pathways, and it coordinately regulated multiple 
developmental processes and stress responses [63]. As 
for AP2-ERF domain containing TFl, its expression pat- 
tern was different from other five members. It showed 
high similarity with DREB26 (AT1G21910.1). In plant, 
RAP2.6, RAP2.6 L, DREB26 and DREB19 exhibited tis- 
sue specific expression and participated developmental 
processes as well as biotic and/or abiotic stress signaling 
[64]. Though previous researches emphasized the func- 
tions of these AP2-ERF TFs on resistance against biotic 
and abiotic stresses, AP2-ERF TFs were also participated 
in plant development such as embryo patterning [65], 
and stamen emergence [66]. 

Additionally, two MYB (R2R3-MYB) transcription fac- 
tors also showed differential expression between QS and 
EG. In plant, MYB TF family was categorized into 3 sub- 
families according to the number of adjacent repeats of 
MYB-domain. Of them, R2R3-MYB subfamily contains 
the largest number of members. Like the AP2-ERF TF 
family proteins, MYB family proteins also function in vari- 
ous plant-specific processes. In Arabidopsis, MYB TFs 
were found as key regulators involved in development, 
metabolism and biotic and abiotic stress responses. 
Among these MYB TFs of Arabidopsis, AtMYB26 is 
involved in determining endothecial cell development 
within the anther and is essential for anther dehiscence 
[67]. AtMYB33 and AtMYB65 redundantly facilitate an- 
ther and pollen development [68]. AtMYB80 regulates 
exine formation and acts downstream of AtMYB35; and 
AtMYB103 is required for tapetal development and 
microsporogenesis, especially for callose dissolution and 
exine formation [69,70]. AtMYB125 positively control 
male germ cell division and commit progenitor germ cells 
to sperm cell differentiation [71,72]. In rice, CSA gene en- 
coding MYB TF functions as a key transcriptional regula- 
tor for sugar partitioning during male reproductive 



development, and the CSA mutant showed reduced levels 
of sugars and starch in floral organs which lead to MS. 

Interestingly, in our results, one MYB TF showed simi- 
lar expression pattern with AP2-ERF TFs that down- 
regulated at BF stage when the anther and pollen grains 
are mature. This MYB TF termed as R2R3-MYB TF was 
closely related to ATMYBR1/ATMYB44 (AT5G67300.1), 
and AtMYB44 was likely to enhance drought and salt 
stress tolerance by suppressing the expression of genes en- 
coding PP2Cs, which was described as negative regulators 
of ABA signaling [73]. Previous report showed that 
AtMYB44 was with changed expression during late em- 
bryogenesis and seed maturation [74]. And notably there 
was a NAC domain protein (JU497421) highly homolo- 
gous with ANAC102 (AT5G63790.1). ANAC102 was an 
important regulator of seed germination and activated a 
seed-specific subset of genes under low-oxygen stress; it 
was also necessary for the viability of Arabidopsis seeds 
following low-oxygen treatment [75]. 

In summary, these results suggested that these AP2- 
ERF TFs and the MYB TF functioned redundantly and 
coordinated with other TFs which involved in the com- 
plex network regulating floral organ development. Fur- 
ther research should emphasize on the isolation of 
proteins interacted with these TFs. 

Conclusion 

An integrative approach combining SSH and microarray 
was employed to explore the transcriptional changes of a 
seedless bud sport mutant of Ponkan mandarin. A num- 
ber of differentially expressed genes were identified. And 
the majority of genes were down- regulated in the mu- 
tant, especially those related to basic metabolic process. 
Metabolism of nutrition and energy might be impaired 
during male gametophyte development of the mutant, 
and TFs and phytohormones might play important regu- 
latory roles during this process. Our research gained 
general information of citrus MS at transcription level 
and could provide some clues for further exploration of 
MS in citrus species. 

Methods 

Accession numbers of sequences and microarray data 

All the sequences generated in the study were deposited 
in GenBank with accession numbers from JU497308 to 
JU497435. Five sequences which are shorter than 200 bp 
longer than 100 bp are attached in Additional file 3. 
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Microarray data and experimental information from 
this study were deposited in the Gene Expression Omni- 
bus (http://www.ncbi.nlm.nih.gov/geo/) under accession 
number GSE38094. 

Plant materials and phenotype analyses 

Two Ponkan mandarin {Citrus reticulata Blanco) culti- 
vars, Qianyang seedless (QS, mutant type) and Egan 
NO.l (EG, a common seedy Ponkan mandarin, wild type) 
were grown in the same orchard of 'Fenghuangshan' citrus 
production area in the city of Dangyang, Hubei province, 
China. These two scion cultivars were seven years old when 
sampling in 2010, with trifoliate orange {Poncirus trifoliata 
L. Raf.) as the rootstock. Flower samples were collected 
from both cultivars in parallel including 4 continuous 
phonologically developmental stages (Figure 7C): squaring 
stage (SF, about 20 DBF), medium bud stage (MF, about 
10 DBF), flowers at full bloom stage (BF) and young 
ovaries of 2-3 days after flowering (OV). All the flowers 
were bagged to prevent cross -pollination, and when 
sampled in the field, all the samples were frozen in liquid 
nitrogen as quickly as possible and then stored at -80°C 
until needed. 

The morphology of mature anthers were investigated 
with fluorescence stereo-microscope (Figure 7A; 7B) 
(Leica MZ FLIII, German) and image was captured with 
a digital camera (Nikon Coolpix, Japan). The pollen 
grain number per anther was counted. In brief, anthers 
from mature flowers were collected and mixed ran- 
domly, each time 40 anthers were dissected and pollen 



grains were suspended in 25 mL sterile water with 4-5 
drops of surfactant (Tween-20, Amresco solon, OH). 
The viability of mature pollen grains were evaluated by 
dying with 1% acetic acid magenta as well as 1% iodine 
potassium iodide (I2-KI2) solution. After staining for 
5 min, pollen grains were observed using BX-61 fluores- 
cence microscope (Olympus, Japan) and Images were 
captured with DP70 CCD digital camera system. At least 
1,000 pollen grains were counted. These experiments 
were repeated three times. The morphology of pollen 
grains was examined by scanning electron microscope 
(SEM) (NTC JSM-6390LV, Japan). For SEM, anthers at 
various developmental stages were pre-fixed with 2.5% 
glutaraldehyde in 0.1 M sodium phosphate buffer (pH 7.2) 
for 24 h, dehydrated twice using a gradient ethanol serial 
(30%-50%-70%-85%-95%-100%), then replaced ethanol 
with isopentyl acetate for 20 min. After that, samples 
were dried with critical-point drying method then 
sputtered coating with gold. Representative images were 
captured. 

RNA extraction and mRNA isolation 

The materials (floral organs) for RNA extraction were 
sampled from at least six independent plants, and mixed 
randomly. Total RNA from flower samples at four stages 
(SF, MF, BF and OV) were extracted with modified Trizol 
method according to [76]. The RNA pellets were washed 
with 75% (V/V) ethanol twice, dissolved in RNase-free 
water and stored at -80°C until use. By mixing equal 
amount of RNA of the four stages, RNA pools from both 



SF MF BF 



Figure 7 Flower organs at different developmental stages and mature anthers. A, B showed the anthers of EG and QS respectively; C 
showed the flower organs including SF, MF, BF stages, the upper row was EG. 
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QS and EG were established in parallel Then mRNA was 
isolated from each of the RNA pools using the Oligotex 
mRNA mini kit (Qiagen, Germany). The quality of RNA 
was determined by Nanodrop 1000 spectrophotometer 
(Thermo Scientific, Wilmington, DE, USA) and 1.2% agar- 
ose gel electrophoresis. 

Suppression subtractive hybridization (SSH) cDNA 
libraries construction and cDNA inserts amplification 

Two micrograms of mRNA was used to synthesize cDNA 
for suppression subtractive hybridization (SSH). The SSH 
was performed with the PCR-select™ cDNA subtraction 
kit (Clontech, Palo Alto, CA, USA) according to the user 
manual. And both forward (the seedless cultivar QS as 
tester and the seedy cultivar EG as driver) and reverse (EG 
as tester while QS as driver) SSH were conducted. For 
cDNA libraries construction, two hybridizations were per- 
formed followed by two rounds of PGR amplifications to 
enrich the desired differentially expressed sequences. 
Then the second PCR-amplified cDNAs were purified and 
ligated into the T/A cloning vector pMD18-T (Takara, 
Japan) overnight at 4°C. Then the ligated products were 
transformed into Electro MAX™ DH5a-E™ cells (Invi- 
trogen, USA) and incubated at 37°C, 160 r/m for 1 h, then 
cultured on SOB-MgCl2 solid media with ampicillin 
(60 [ig ml'^) to generate the primary cDNA libraries. The 
transformed white bacteria were randomly picked and 
grown on 384-well plates containing Luria Broth (LB) li- 
quid media with ampicillin (100 [ig ml"^) at 37°C overnight 
(about 16 h). Glycerol (Amresco, USA) (4.4% final) was 
added for storage at -80°C. 

A total of 8,000 cDNA clones were randomly picked 
from forward and reverse SSH libraries and used as for 
subsequent PGR templates. Each PGR was performed in 
a 100 [A reaction mixture using nested primers of SSH 
according to [77]. The PGR products were precipitated 
with equal amount of isopropyl alcohol and washed with 
75% (V/V) ethanol, then re-suspended in 40 [A sterile 
water. The yield and quality of the PGR products were 
determined by Nanodrop 1000 spectrophotometer 
(Thermo Scientific, Wilmington, DE, USA), and then run 
on 1.2% agarose gel and examined by Bio-Rad UV spec- 
troscopy (Bio-Rad Laboratories, Washington, DG, USA) 
to confirm single clone (Additional file 4: Figure S3). Fi- 
nally the validated PGR products were stored at -80°G 
for custom microarray. 

Microarray slides fabrication and preparation of 
fluorescent dye-labelled cDNA 

About 40 microlitre of PGR products were re-precipitated 
by adding 100 [A of anhydrous ethanol and were dissolved 
in EasyArray^^ spotting solution (GapitalBio Gorp, Ghina) 
at a final concentration of 0.1-0.5 [ig [A'^ and then printed 
on amino-silaned glass slides with a SmartArrayer^^ 



microarrayer (GapitalBio Gorp). Each clone was printed 
triplicate. The particular procedures for microarray fabri- 
cation were conducted according to [37]. 

The relative gene expression profiles of QS at four de- 
velopmental stages (SF, MF, BF and OV) compared with 
the corresponding four stages of EG were investigated 
by microarray analysis. For each stage, three sets of total 
RNA samples were extracted independently, and then 
RNA pool was constructed by mixing aliquot of RNA 
from the three sets of RNA samples. An aliquot of 5 [ig 
total RNA from the RNA pool was used to produce 
Gy5/Gy3-labelled cDNA employing an RNA amplifica- 
tion combined with Klenow enzyme labeling strategy 
according to the protocol by [78]. Gy5/Gy3 -labelled 
cDNA was hybridized with the microarray at 42°G over- 
night. Hybridization was performed in duplicate by dye 
swap (Gy5-labelled cDNA of QS versus Gy3-labelled 
cDNA of EG, and Gy5-labelled cDNA of EG versus Gy3- 
labelled cDNA of QS). And then the arrays were washed 
with 0.2% SDS, 2 x SSG at 42°G for 5 min, and 0.2% SSG 
for 5 min at room temperature. 

Microarray data analysis and EST sequence analysis 

Arrays were scanned with a confocal laser scanner, 
LuxScan™-scanner (GapitalBio Gorp.) and the resulting 
images were analyzed with LuxScan^^ 3.0 software 
(GapitalBio Gorp.). cDNA spots were screened and iden- 
tified with the methods described by [77]. A spatial and 
intensity-dependent (LOWESS) normalization method 
was employed and normalized ratio data were then log 
transformed [79]. Differentially expressed genes were 
identified using a t-test, and multiple test corrections were 
performed using FDR. Genes with FDR <0.05 and a fold 
change >2 were identified as differentially expressed 
genes. 

All the clones differentially expressed in at least one of 
the four stages were subjected to single-pass sequence 
using standard high throughput sequencing by BGI- 
Wuhan, Ghina. All sequences were edited to omit vec- 
tors and low quality segments at 5' and 3' ends, then re- 
moval of sequences shorter than 100 bp with SeqGlean 
software. Sequence reads were assembled by GAP3 pro- 
gram [80] with default parameters. Then all the unigenes 
were annotated using BLASTx with a cut-off value of 
1.0 X e"^ by searching the UniProt database (http://www. 
ebi.ac.uk/uniprot/). GO-KEGG-EG annotation was per- 
formed based on Annot8r platform [81]. Hierarchical 
clustering of transcript accumulation was performed 
with Gluster software (version 3.0) [82]. 

Quantitative real-time PGR verification and candidate TFs 
analysis 

Total RNA was extracted from QS and EG collected at 
four different developmental stages with the Trizol 
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methods mentioned above. Primer pairs were designed 
with the Primer Express software (Applied Biosystems, 
Foster City, CA, USA). Primer sequences of 11 candidate 
genes for verification were provided in Additional file 5: 
Table SI, and primer sequences of 10 TPs were provided 
in Additional file 6: Table S2. Single strand cDNA was 
synthesized with the prescription of the Revert Aid TM 
first strand cDNA synthesis Kit (Fermentas, Life Science, 
EU). Then each cDNA sample was pre-amplified using 
the citrus house-keeping gene |3-actin and normalized 
for subsequent real-time quantitative PGR (qRT-PCR). 
The PGR program differed in terms of the annealing 
temperature of each primer pair and the length of the 
predicted PGR products. The qRT-PGR was per- 
formed using the ABI 7500 Real Time System (PE 
Applied Biosystems, Foster Gity, GA, USA) with the 
method as described by [83]. And relative transcript 
change was analyzed by 2'^^^^^\ 

Additional files 



Additional file 1: Figure SI. Comparison of expression patterns of 4 
ESTs assembled the some contig encoding Methionine synthase. Y-axis 
represents the average Cy5 (Cy3) to Cy3 (Cy5) ratio in array hybridization. 
X-axis represents the four developmental stages. 

Additional file 2: Figure S2. The correlation of gene expression ratios 
between cDNA microarray and qRT-PCR. Data were from 1 1 probe sets at 
four developmental stages. The gene expression ratios based on cDNA 
microarray were in log2 transformed. 

Additional file 3: Sequences shorter than 200 bp but longer than 
100 bp. 

Additional file 4: Figure S3. The purified PCR products for microarray 
probe. 100 bp molecular ladders were used. 

Additional file 5: Table SI. qRT-PCR primers for 1 1 candidate genes 
and citrus actin gene. 

Additional file 6: Table S2. qRT-PCR primers for 10 transcription factors 
aFs). 
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